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ABSTRACT 

We present in this article an analysis of some of the properties of the density field 
realized in numerical simulations for power-law initial power-spectra in the case of a 
critical density universe. We study the non-linear regime, which is the most difficult 
to handle analytically, and we compare our numerical results with the predictions of 
a specific hierarchical clustering scaling model that have been made recently, focusing 
specifically on its much wider range of applicability, which is one of its main advantages 
over the standard Press-Schechter approximation. We first check that the two-point 
correlation functions, measured both from counts in cells and neighbour counts, agree 
with the known analytically exact scaling requirement (i.e. only depend on ct^) and 
we also find the stable-clustering hypothesis to hold. Next we show that the statistics 
of the counts in cells obey the scaling law predicted by the above scaling model. 

Then we turn to mass functions of overdense and underdense regions, that we 
obtain numerically from "spherical overdensity" and "friends-of-friends" algorithms. 
We first consider the mass function of "just-collapsed" objects defined by a density 
threshold A = 177 and we note as was found by previous studies that the usual Press- 
Schechter prescription agrees reasonably well with the simulations (although there are 
some discrepancies). On the other hand, the numerical results are also consistent with 
the predictions of the scaling model. Then, we consider more general mass functions 
(needed to describe for instance galaxies or Lyman-a absorbers) defined by different 
density thresholds, which can even be negative. The scaling model is especially suited 
to account for such cases, which are out of reach of the Press-Schechter approach, and 
it still shows reasonably good agreement with the numerical results. Finally, we show 
that mass functions defined by a condition on the radius of the objects also satisfy the 
theoretical scaling predictions. 

Thus, we find that the scaling model provides a reasonable description of the 
density field in the highly non-linear regime, for the cosmologies we have considered, 
both for the counts in cells statistics and the mass functions. The advantages of this 
approach are that it clarifies the links between several statistical tools and it allows 
one to study many different classes of objects, for any density threshold, provided one 
is in the fully non-linear regime. 
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1 INTRODUCTION 

In the standard cosmological scenario large-scale structures 
in the universe arise from the amplification through gravi- 
tational instability of small primordial density fluctuations. 
These initial perturbations are likely to be gaussian (as in 
most inflationary models) and are characterized by their 



power-spectrum. Within the hierarchical clustering scenario, 
the amplitude of these fluctuations increases towards the 
smaller scales (as in the CDM model: Peebles 1982; Davis 
et aI.1985). Small scales collapse first to form bound ob- 
jects which merge later to build increasingly massive halos 
as larger scales become non-linear. These mass condensa- 
tions correspond to the various astrophysical objects one 
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can observe in the universe, from Lyman-a clouds to galax- 
ies and clusters. As a consequence, to understand the for- 
mation of these objects it is important to obtain a precise 
description of the evolution of the density field under the 
action of gravity. Note that to model specific astrophysical 
objects one usually needs to add to this description of dark 
matter halos non-gravitational effects (such as, for instance, 
radiative processes to get the evolution of baryons). This is 
beyond the scope of the present paper that considers specif- 
ically the multiplicity of dark halos. 

However, even for the problem of the evolution of the 
matter distribution under the sole action of gravity theo- 
retical results are very scarce. Linear theory allows one to 
describe the density field on large scales, while a few approx- 
imations (e.g. see Bernardeau 1996) try to handle the early 
non-linear evolution, but the highly non-linear regime has 
proved very difficult to model. As a consequence numerical 
simulations have so far been the main tool to describe this 
latter stage. In this article, we compare the results of simu- 
lations for power-law initial spectra in the case of a critical 
density universe with the scaling model presented in previ- 
ous publications (Valageas & Schaefi'er 1997, hereafter VS; 
Bernardeau and Schaeffer 1991; Balian & Schaeffer 1989a ). 
This latter description of the density field is based on the 
assumption that the many-body correlation functions sat- 
isfy specific scaling laws (obtained from the stable-clustering 
ansatz) in the highly non-linear regime. In this case one ob- 
tains a very powerful model for the density field, which can 
be used to obtain the counts in cells statistics as well as mass 
functions of overdense and underdense regions at difi'erent 
density thresholds. Note that in contrast the Press-Schechter 
(PS) mass function (Press & Schechter 1974), which is cer- 
tainly the most popular tool to get some information on the 
characteristics of the non-linear density field, only deals with 
"just-collapsed" objects. 

The article is organized as follows. In Sect.H we describe 
the numerical simulations we analyse. Then we present our 
results for the two-point correlation function and the counts 
in cells statistics in SectM We compare the numerical mass 
function of "just-collapsed" objects to the usual PS prescrip- 
tion and to the scaling model in Sect.M. We also consider 
more general mass functions, beyond the reach of the PS 
approach, that are defined by various density thresholds. 
Finally we present our results for the limiting case of mass 
functions of objects defined by a constant radius constraint. 
These studies originate from astrophysical descriptions of 
galax;ies (Valageas & Schaeffer 1999; Valageas & Silk 1999) 
and Lyman-Q clouds (Valageas et al.l999) where one is nat- 
urally led to introduce such generalized mass functions. 



2 SIMULATIONS 

We shall consider a critical density universe, fl = 1, with an 
initial power-spectrum P{k) which is a power-law: P{k) oc 
fc". We study the cases n = —2, n — —1 and n = 0. As 
usual we shall define a^{R,a) to be the amplitude of the 
density fiuctuations in cells of physical radius R at time t 
(scale-factor a) given by linear theory. Thus we have: 



2 2 -(n+3) (n+5) r>-("+3) 

where r — R/a is a comoving scale. 



(1) 



The N-body simulations were performed using the 
AP^M code of Couchman (1991). The nominal box size 
was L — 256 Mpc/h (though this could be rescaled to any 
value, since the initial conditions are scale-free). The simu- 
lations all used Np = 128'^ ~ 2 x 10^ particles, with a force- 
softening parameter (constant in comoving coordinates) of 
0.1 times the mean interparticle separation L/Np , that is 
over a comoving radius of 0.2 Mpc/h. The expansion factor 
a was normalized so that (j{R, a) = 1 for R — 8 Mpc/h 
at a = 1, according to linear theory. The initial positions 
and velocities of the particles were given by displacing the 
particles from a cubical grid using the Zeldovich approxi- 
mation and the linear theory power spectrum. The starting 
time was chosen small enough so that the density fluctua- 
tions on the scale of the particle grid were still close to the 
linear regime. Specifically, the linear power-spectrum ampli- 
tude at the Nyquist frequency of the particle grid was chosen 
to be A^ times the white noise value, with A = 1 for n = 
and n = —1, and A — 0.4 for n = —2, corresponding to 
expansion factors at = 0.0611,0.15,0.194 for n — 0,-1,-2 
respectively. For n — 0, the simulation was evolved with a 
constant timestep Ap = 2.36 x 10^'^ in the variable p — a^''^, 
while for n — —1 and —2 a constant timestep in a was used, 
with Aa = 3.4 x 10"^, 1.5 x 10~^ respectively. The simu- 
lations were evolved up to expansion factors a/ = 8, 4, 2.67 
for n = 0, —1, —2 respectively. In the numerical analysis we 
shall use the output times in the range 1 < a < a/ when 
some non-linear structures have already formed on scales 
larger than the smoothening length. 

For this cosmology and initial power spectra, the real 
clustering evolution should be self-similar, when scaled to a 
radius R,{a) oc a("+5)/("+3) such that a{R,,a) = 1. This is 
an exact analytical result that holds independently of the 
validity of the Balian & Schaeffer (1989a ) scaling pre- 
dictions (see Sect.b|) that we aim at testing in this paper. 
The evolution in the N-body simulation will depart from ex- 
act self-similarity because of numerical effects, in particular, 
particle discreteness, force softening, the finite box size, and 
the absence of initial fiuctuations smaller than the Nyquist 
wavelength or larger than the box size. 



3 COUNTS-IN-CELLS 

A convenient way to describe the density field obtained in 
a numerical simulation, or realized in the actual universe, is 
to consider the counts in cells. Thus, we define the proba- 
bility distribution Pii[N\t) to be the probability to have A'^ 
particles in a spherical cell of radius i? at a given time t. In 
the following we shall usually denote PR,{N;t) as P{N) to 
simplify the notation. This is a well-defined quantity which 
provides a very good description of the density field and is 
a convenient tool for a theoretical analysis. In contrast, the 
multiplicity functions used to recover the counts of virial- 
ized halos or astrophysical objects may be defined in various 
ways and are somewhat more difficult to handle analytically. 
However, in certain regimes their properties can be obtained 
from the characteristics of the counts in cells. Hence we shall 
first consider the statistics of P{N) in the next sections. 
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3.1 Non-linear scaling model 

Since we shall mainly compare our results with the scaling 
model described in Balian & Schaeffer (1989a ) we recall here 
their predictions while introducing our notation. This model 
is based on the assumption that the many-body correlation 



functions ^p(ri, ... 
^p{Xri,...,Xrp-a) 



Vp ; a) satisfy the scaling law : 



3(p-l) ^-7(P-1) 



^p(ri 



(2) 



where a{t) is the scale-factor and 7 is the slope of the two- 
point correlation function (which we note ^). For an initial 
power-spectrum P{k) which is a power-law, P{k) ex k", we 
have: 

-^ ,3, 

The relations (bl) and (pi) are derived from the stable- 
clustering assumption (Davis & Peebles 1977; Peebles 1980) 
and are expected to describe the highly non-linear regime 
^ ^ 1. To obtain the statistics of P{N) in cells of radius R, 
volume V, it is convenient to define the quantities: 

S. = 4^ with f^ = / ^%^ Uru-:rp). (4) 



e' 



Vp 



Then, 

ptN) = — L h(x) with x= ^^ -!-4^ (5) 

in the regime where 

iV>iV„ and iV^ > iV„. 

Here we defined 



(6) 



Nc = N^, N, 



N 



-c^/(l-^)' 



N = nV, 



N 
1 + A = = (7) 

iV 



where n is the average number density in the simulation. It is 
also required that the continuous limit A*' ^ 1 is reached in 
the relevant sampling of the simulation. The two conditions 
(a) insure that i) the counts are relevant to mass condensa- 
tions and not to "voids" and ii) one is in the fully non-linear 
regime where ^ is large so that the non-linear scaling repre- 
sented by h{x) can develop. The scaling (H) implies that 



p > 1 : Sp = / x^hix) dx , S'l = S2 = 1. 
Jo 

From very general considerations one expects that 
s < 1 : h{x) 



a{l-uj) ^_2 

X 



T(lo) 



(8) 



(9) 



X ^ 1 : nix) ^ as X e ' 



where Xs is a constant, typically of order 10. Thus, P{N) 
shows a power-law behaviour in the range N^ <gi N <gi Nc 
and an exponential cutoff for N ^ Nc- 

If the scaling laws (S) apply, then the ratios Sp and the 
function h{x) are independent of scale and time. This means 
that once h{x) is given (e.g. from measures performed at a 
certain scale and time) one only needs to know the evolu- 
tion of the two-point correlation function ^ to be able to 
construct the whole statistics of P{N) at any scale and time 
in the highly non-linear regime. Since $, obeys (^j , its evolu- 
tion is known in the non-linear regime once its normalization 
is measured for one scale and time. 



Table 1. Normalization parameters for the two-point correlation 
functions. The non-linear exponent 7 is given by (pi). The columns 
(J) and (C) correspond to the numerical results obtained by Jain 
et al.(1995) and Colombi et al.(1996), see text. Note that the ratio 
a/a is predicted by the stable-clustering ansatz, see (jig). 



n 


7 


a/a 


a 


a (J) 


"(C) 


a 


a (J) 


"(C) 





1.8 


1.10 


1.71 


1.43 


1.28 


1.55 


1.30 


1.16 


-1 


1.5 


0.95 


1.45 


1.18 


1.33 


1.52 


1.24 


1.40 


-2 


1 


0.96 


1.25 


1.09 


1.25 


1.30 


1.13 


1.30 



3.2 The two-point correlation function 

3.2.1 Fluctuations in a cell 

As we described in the previous section, in order to test 
the scaling model in the domain A'' S> Nv we first need the 
evolution of the two-point correlation function ^. Moreover, 
^ presents a strong interest in itself since it gives a measure 
of the amplitude of the density fluctuations and it provides 
a first check of the stable-clustering assumption. Indeed, if 
the latter holds the slope 7 of ^{R) must be given by (H). 

To get ^{R) we simply count the number of particles 
enclosed in each of 300"^ spheres of radius R set on a grid 
and we measure (A'') and {N'^), where {) denotes an average 
over all trials. Then we use the relation: 



c 






1 



1 

JN) 



(10) 



As argued by VS, following the ideas of Hamilton et 
al.(1991) from a Lagrangian point of view where one follows 
the evolution of matter elements, one expects ^{R, a) to be 
closely related to the linear correlation function a^{RL,a) 
evaluated at a different scale Rl but at the same time: 



^iR,a) = F[a\RL,a)] 
Rl = [l+l{R,a)] R' 



(11) 



In fact, since our initial conditions are scale-invariant the 
relation (hll) is exact: it is then just a rewriting of the scal- 
ing solution of Peebles (1980) that holds in this case. Its 
main interest comes from its physical interpretation which 
suggests that the F{a'^) obtained for various n should show 
a similar behaviour. Of course in the linear regime ^ <C 1 we 
must recover ^ ~ ct^ while in the highly non-linear regime 
if clustering is stable we have ^{R,a) oc a^. This gives the 
asymptotic behaviour of the function F{a'^): 



o- < 1 : F{a^) = a^ 



cj:» 1 



n-') = m -'- 



(12) 



where q ~ 1. As a consequence, we plot the values of £,{R) 
obtained from the simulations through (y_0|) as a function of 
a^{RL) rather than R. Thus, at a given scale R we measure 
£,{R) and we know (J^{R) from which we derive (j^{Rl) = 



0-^(7?) {Rl/R) 



-(n + 3) 



a^{R) (1 + C)"'"+^'^^- The results 



are shown in Fig.g] for various n. 

Since the initial conditions of the simulation are scale- 
invariant {P{k) is a power-law) all curves should superpose 
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Figure 1. The two-point correlation function ^{R) as a func- 
tion of a^{Ri,) for the power-spectra n = 0,-1 and —2. The 
sohd hne is the asymptotic behaviour ( |l2[ ) for large ^, in the 
stable clustering regime, with normalization a from Tab.l. The 
crosses are numerical values obtained from counts-in-cells, while 
the squares are the estimates of § provided by measures of ^ from 
neighbour counts, see main text and (|l5[). Different shades of grey 
correspond to different comoving scales (0.2, 0.5, 1, 2, 4, 8 and 
16 Mpc). The filled (resp. open) triangles show the results for 
0.2 Mpc from counts-in-cells (resp. neighbour counts). The larger 
scales (8 and 16 Mpc) saturate at too low a value of 5, reflecting 
the finite size of the sample. The 0.2 Mpc scale presents devia- 
tions from the scaling in the highly non-linear regime (all curves 
should exactly superpose for our power-law initial conditions), 
with a lack of power due to the softening-over a 0.2 Mpc radius- 
of the gravitational interaction. 



as long as numerical effects (finite box size, softening length 
and resolution scale) are small. We can check that this is 
indeed the case. Moreover, in the highly non-linear regime 
we recover the asymptotic behaviour given by ([12) which is 
shown by the solid line. The normalization parameter a is 
displayed in Tab.l for the three power-spectra. However, we 
can note that for larger scales the two-point correlation func- 
tion "saturates" to a value lower than its asymptotic limit. 
This problem increases for smaller n and is due to numerical 
effects since as we explained above all curves should super- 
pose. The dependence on scale and on n of this effect shows 
that it is produced by the lack of power in the simulation 
at large scales. Indeed, the actual initial power-spectrum 
P{k) is a power-law only over a limited range of k and the 
influence of the numerical cutoff at kmin, unduely supress- 
ing some power, becomes more important as one considers 
larger scales (hence wavenumbers closer to kmin = 2tt/L) 
and lower n (which increases the contribution of small k). 
Similarly, there is a high frequency cut-off. The correlation 
function ^ starts to deviate from the exact scaling imposed 
by our power- law initial conditions when the sampling of the 
numerical output is done at the 0.2 Mpc comoving scale and 
below. This is very distinctly seen in Fig.hl for the larger val- 
ues of the correlation function. This is due to the softening 
of the interaction, which for all the samples we use is done 
over a radius of 0.2 Mpc (comoving). Again, the result is a 
lack of power, due to this softening. Obviously, the deficit 
is larger for n = than for n = —2, the former case, with 
more small-scale power, being more sensitive to this effect. 
Note that the Nyquist frequency and particle discrete- 
ness play a minor role if any at all since we only consider the 
cases where the number of points in the simulation is large 
enough. 



3.2.2 Counts of neighbours 

Although the average ^ defined by (W) appears naturally 
when one considers counts in cells many authors studied 
the quantity ^ given by: 



m - ^ 



^(r) r dr 



(13) 



which is relevant for the counts of neighbours. Indeed, the 
mean number of neighbours located within a radius R from 
a given particle is: 



(iV„,) = {N) (l + |(i?)) 



(14) 



For a power-law power-spectrum ^ depends on its lin- 
ear equivalent a^ in a fashion similar to (hll) with a new 
function F which has the same asymptotic behaviour as 
([12) but with a slightly different normalization a (Hamil- 
ton et al.l991; Peebles 1980; see also Padmanabhan 1996). 
We measure ^ in the simulation using (H): we count the 
mean number of neighbours of particles in the box which 
gives < Nnb > while the total number of points in the sim- 
ulation provides < N >. The result is shown in Fig.0. We 
indeed recover the behaviour predicted by (|l3). This con- 
firms the previous results of Fig.|l| since both ^ and ^ are 
measures of ^. This shows that in the non-linear regime the 
stable-clustering ansatz holds for ^ which obeys the scaling- 
law (0). 
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Figure 2. The two-point correlation function §(-R) (obtained 
from neighbour counts) as a function of a^{R]^) for the power- 
spectra n = 0,n = — 1 and n = — 2. The soUd hne is the asymp- 
totic behaviour (tl3) for large ^, with the normalization parameter 
a from Tab.l. The triangles correspond to the comoving scale 0.2 
Mpc. See FigJll and the text for comments. 

Moreover, the measure of ^ allows us to get a second 
estimate of ^ in the non-linear regime. Indeed, in the domain 
where ^ ^ 1 the two-point correlation function is a power- 
law with a slope 7 given by (H) as shown by Fig. hi and Fig.bl 
In this case, one can show (Peebles 1980) that: 



|(i?) = /? C(i?) with /? = (1 - 7/4) (1 - 7/6) 2^ (15) 
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Thus, we can obtain ^{R) from the measure of £,{R)- These 
values are shown by squares in Fig.ty. We can check that 
they agree with the previous calculations for £, based on the 
counts in cells. Note that we only show in FigJll the values 
of ^ obtained from ^, through (\laj), which are in the highly 
non-linear regime ^ 3> 1 where the two-point correlation 

applies. Note that in 
, the lack of power at 



function is a power-law so that (tl! 
this case (compare Fig.nl with Fig! 
the larger scales is much less pronounced: the count of neigh- 
bours already focuses on the dense regions and the statistics 
of the poorly represented events at large scale is improved. 
Indeed, to measure ^ one considers cells which are centered 
on particles, hence this statistical tool follows the evolution 
of gravitational clustering as it automatically probes more 
closely denser regions, while to get ^ the center of the cells 
is set at random in the box, so that in the highly non-linear 
regime most cells are within voids. Thus at least part of 
the deviations seen at these scales are due to the statistical 
extraction of the information. 

We present in Tab.l the values of the normalization 
parameters a and a obtained from our simulation, as well 
as from Jain et al.(1995) and Colombi et al.(1996). For the 
former we calculate a from their value of a while for the 
latter we obtain a from a. Indeed, from (hsl) we get: 



a (i^l^ /3-i/^ 



(16) 



where /3 (resp. /3l) is given by (113) with 7 — 3(3-|-n)/(5-|-n) 
(resp. 7 = (3-|-n)), corresponding to the non-linear (resp. lin- 
ear) regime. Although all numerical simulations agree with 
the stable-clustering ansatz: (|3|), ([ll|) and (|l2|), the normal- 
ization of ^ in the non-linear regime varies up to a factor 2 
(note that it scales as the cube of a). Thus, there is still some 
inaccuracy in the numerical values of ^. Moreover, one may 
expect the higher-order correlation functions (hence the pa- 
rameters Sp) to bear at least similar uncertainties. The dis- 
crepancies between the various estimates of a may be due 
to the effects of finite volume and particle number: in par- 
ticular the initial power-spectrum is not a power-law over 
an infinite range of scales (there are an upper and a lower 
cutoff) and one should average ^ over many realizations of 
the initial gaussian field. However, the fact that curves ob- 
tained from different scales superpose in Fighl and Figj2| and 
that our measures of ^ are consistent with ^ show that we 
can reasonably rely on our results. Except in the extreme 
cases where the finiteness of the sample (at large scales) or 
the smoothing of the force (at small scales) become impor- 
tant, these results are consistent with the idea that the stable 
clustering regime is reached (see however Padmanabhan et 
al.l996 for an alternative view). 



3.2.3 Analytical fit to the fluctuations in a cell 

We show in Fig.H the relation cr^{R) «-» £,{R). This cor- 
responds to an Eulerian point of view as opposed to the 
Lagrangian point of view displayed in Fig.nl, since we now 
consider the values of the linear and non-linear correlation 
functions taken at the same spatial scale (and not at the 
same "mass scale"). We can note that the curves we ob- 
tain are quite smooth when displayed in this diagram and 
they do not show a sharp "step-profile" as was the case in 
Fig.hl and Fig.H. The solid lines are analytical fits to the 
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Figure 3. The two-point correlation function i;{R) as a function 
of o-^{R) (hence taken at the same spatial scale). The crosses are 
numerical estimates from counts in cells in the simulations while 
the solid lines are the analytical fits \vA . See Fig Jll and text for 
comments. 

numerical results which we shall need below (SectJ^ to get 
the mass functions within the framework of the scaling ap- 
proach. Thus we use the approximations: 

' n = o,i: aR) ^ [i^')"' + iiirT'^'^ 

with n = : a = 0.7 ; n = -l:a = 4 (17) 

n = -2: aR) ^ [<^' + i'^V5)%] I [1 + (aV5)'] 



where i^{R) = (|£ rj{R)fl'-'^+''\ Of course, these fits are 
built so as to be consistent with (llll) and ([12|). These re- 
lations define implicitly the functions F(a ). They will be 
used instead of the correlation function E, actually measured 
in the simulation (hence they will not be a source of un- 
wanted deviations from self-similarity) when we determine 
the scaling properties of the multiplicity functions (Sect.l4|). 



3.3 Counts in cells statistics 

By counting the number of particles enclosed in each of 300"^ 
spheres set on a grid, as we did to evaluate ^, we can also 
obtain the statistics of the counts in cells. This allows us 
to compare the predictions (H) of the scaling model to our 
numerical results for P{N). Note that we obtain simultane- 
ously the value of ^ in the simulation, at the scale of interest, 
which we use in (pi) . We display the results in Fig.W as a func- 
tion of x: from A^ and P{N) obtained by these counts in the 
simulation we define; 

t2 






and X h{x) 






P{N) 



(18) 



If the scaling-laws (Q) hold the quantity h{x) used above 
must be the scaling function defined in (H) . Then, all curves 
obtained for different sizes and times should superpose. Note 
that in our case where the initial power-spectrum is a power- 
law all curves characterized by the same ^ (or equivalently 
by the same cr^) must coincide. Thus the scaling-laws (bh 
merely imply in addition that curves measured for different ^ 
should also superpose, when shown in a diagram x ^^ x^h{x) 
as defined by wM- 

We can see in FigMl that the various curves indeed su- 
perpose, for all three power-spectra. Note that in this fig- 
ure we only show the parts of P{N) which should scale as 
predicted by (H): we only display points which satisfy the 
requirements ^ > 100, N > 4,N > 4iV„ and Nc > 4. The 
solid curves are analytical fits to the data points. We use the 
functional form introduced by Bouchet et al.(1991): 



h{x) = 



a(l 



exp(— x/a;s) 



(19) 



T{uj) {1 + bxy 

The values of the parameters a,uj,b,c and Xs are given in 
Tab. 2. Note that the functions h{x) must also obey the con- 
straints: 



X h{x) dx = 



x h{x) 



dx = 1 



(20) 



from the general relation (pi). 

The dashed curve in the figures is the function h(x) 
which one obtains (see VS) by assuming that the non-linear 
evolution of initial density fluctuations is given by the Press- 
Schechter approximation in the case where the correlation 
function has reached the stable clustering regime (which is 
indeed the case in the examples we have considered in FigM). 
Thus, we obtain what we shall refer to as the "PSs" estimate 
("s" for stable) 



hpss{x) 



2 5 + n ^(„+5)/6-2 ^^p 
TV 6a 



j.(5-l-n)/3 

2a2 



(21) 



the result being multiplied by the standard "PS factor" of 
two. For the stable clustering regime, a; < (1 + A)/200, this 
expression is identical to the usual Press- Schechter formula 
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Table 2. Parameters for the scaling function h{x). 



log(x) 

Figure 4. The probability distribution P{N) for counts in cells. 
We display the quantity N^/N P{N) = x^h{x) as a function 
of X = N/Nc- The solid curves are the analytical fits given by 
( |l9[ ) and Table 2 while the dashed lines show the "PSs" scaling 
function (bll) obtained using the Press-Schechter approach in the 
regime where the correlation function has reached stable cluster- 
ing. Different shades of grey correspond to different values of §, 
ranging from 100 up to the highest value shown in Fig.[ll. Only 
the scales (from 0.5 to 4 Mpc) where the correlation function 
is seen (Figs.hl, H, H) to obey the scaling required by the exact 
self-similarity are considered here and in the following. 



n 


a 


UJ 


Xs 


b 


c 





1.50 


0.45 


4 


1 


0.6 


-1 


1.48 


0.4 


8 


3 


0.6 


-2 


1.71 


0.3 


13 


5 


0.6 



(normally used for A ~ 177, the usual density contrast of 
just-virialized halos). If the above condition on x is not ful- 
filled, only minor differences between the two expressions 
appear however, see FigW below. We can see in Fig.W that 
the function hpss{x) obtained in this way does not agree 
with the numerical results, the true distribution being much 
broader as argued by VS (the general trend, however, larger 
ui and sharper cutoff for higher n is correct) . This means ei- 
ther that one cannot follow the evolution of individual mat- 
ter elements with using only the spherical model, because 
there is some exchange of matter between neighbouring den- 
sity fluctuations through mergers and disruptions and the 
corrections to the spherical dynamics are not negligible, or 
that after virialization the density fluctuations still undergo 
non-negligible evolution. The stable-clustering ansatz on the 
other hand is seen from Figjlj and Figjj to be a valid approx- 
imation in a statistical sense. 

The scaling of the count-in-cells has thus been checked 
over an unprecedented range (more than 3 decades in the 
scaling parameter x). Note that we check here not only the 
scaling of the counts as a function of the cell size R but also 
as a function of time, that is the self-similarity -in principle 
exact- of the computation results. 

Although all curves at different scales superpose, build- 
ing indeed the predicted universal h{x) curve, some devia- 
tions are seen for small values of the parameter x. Clearly, 
these small values of x (log(a::) = —2.5 !), never reached in 
earlier simulations, are at the limit of the possibilities of 
the present calculations, but the deviations appear to be 
systematic. For larger ^ the function h{x) appears to be 
slightly flatter although it still obeys the constraints (EOl). 
Whether this effect is the sign of a real deviation of h{x) 
from the predicted scaling, or is due to some bias in the 
simulation, warrants investigation. To this purpose, we can 
consider the same simulation results in more detail. In the 
continuous limit (A'^ 3> 1), the expression N^ / N P{N) that 
in general depends on the three variables N , the cell size 
R and the expansion parameter a, can always be written 
as N'^/N P{N) — x^h{x,^,R), the latter expression repre- 
senting a simple change of the original variables into x, ^ 
and R. The exact self-similarity due to our power-law ini- 
tial conditions implies that N /N P{N), as well as 5, are 



a function of the ratio R/i 



,{n + 5)/(n + 3) 



only (i.e. they only 



depend on a) and thus that N^/N P{N) = x^h{x,£,) is 
independent of R since ^ carries all the dependence on R 
and a. The scaling prediction (due to the stable-clustering 
model) in addition states that in the limit of large ^, when 
the stable clustering regime is reached, N"^ /N P{N) has a 
limit for fixed x that is independent of ^. By definition of 
h{x) this limit is iV^/A P(A) = x^h[x). The scatter seen 
in Fig.W must have its origin in the violation of one of the 
two above conditions that are needed to get the seeked scal- 
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log(x) 

Figure 5. Upper panel: the probability distribution P{N) for 
counts in cells as in Fig.W for n = with the constraint 800 < 
^ < 1600. We show the curves obtained for the three comoving 
scales R = 0.5 Mpc (triangles), R = 1 Mpc (squares) and R = 2 
Mpc (circles). The solid curve is the analytical fit given by ( |l9[ ) 
and Tab. 2 as in Fig.kl Lower panel: the curves for P{N) obtained 
at fixed comoving R by letting § (i.e. time) vary. We show the 
cases _R = 1 Mpc (triangles) and R = 4 Mpc (squares). Different 
shades of grey correspond to different values of ^. 



ing function. As a first check, we can consider tlie variations 
witli R at fixed ^. The value of the latter ranging from 100 
to 6400, typical intermediate values are 800 < £, < 1600. 
In this range, values of 7? = 0.5, 1 and 2 Mpc are available. 
We can see in FigH (upper panel) that for these given val- 
ues of ^, although the curves look reasonably close, there is 
a distinct drift at small x, the curves for the smaller val- 
ues of R lying systematically above the ones for the lower 
values. Also, the deviations get larger the smaller the value 
of X. Since we let R vary for a fixed ^ this effect represents 
deviations from the "exact" self-similarity (implied by scale- 
invariant initial conditions) due to the unavoidable errors in 
the simulation. To estimate the size of these deviations, we 
show in FigJa (lower panel) the probability distribution ob- 
tained for all possible values of ^ keeping R fixed at two 
values (1 and 4 Mpc). For a given R, when ^ is varied all 
curves superpose. But when R is changed, we get a differ- 
ent curve. This shows again that N^ /N P{N) depends on 
R and not on ^, with a deviation that reaches a factor of 
2 at small x {x < 3 W~ ), of the order of magnitude of 
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Figure 6. Comparison of the scaling function x^ h(x) obtained 
from our simulations with the estimates given by Colombi et 
al.(1997) (dashed line) and Munshi et al.(1999) (dot-dashed line). 
For the latter two cases we only display x"^ h(x) in the range where 
it was actually tested against numerical data. 

the scatter seen in Fig.W. So we can attribute the observed 
scatter to deviations from the "exact" self-similarity (due to 
a power-law initial power-spectrum) , hence to numerical in- 
accuracies and not to a violation of the non-linear scaling 
prediction that we aim to test here. 

We compare in Fig.g the scaling functions h{x) we get 
in our simulations with the results obtained by Colombi 
et al.(1997) and Munshi et al.(1999). For these two latter 
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cases we only display h{x) over the ranges where there were 
actually data points from the simulations (for our scaling 
function the range tested against numerical results can be 
directly seen from FigMh. We can check in the figure that 
all of the scaling functions agree fairly well with each other 
(within the numerical uncertainties which can be estimated 
from the dispersion shown in Fig.0) . Thus, although there is 
still some dispersion between the various estimates of h{x), 
which could be expected from the discrepancies which al- 
ready appeared for ^, the scaling function h{x) is rather 
well constrained in the range where it has actually been 
tested, to an accuracy which is certainly sufficient for prac- 
tical purposes. However, it should be noted that the range in 
X available to constrain h{x) is still too small to determine 
accurately the parameters of the fit (|l9|) for low n, especially 
Xs (there is some degeneracy between the exponential cut- 
off and the exponent lus of the power- law pre- factor). Thus, 
in the case n = —2 the values of Xs obtained by Colombi 
et al.(1997) and Munshi et al.(1999) diflter by a factor 2, al- 
though both scaling functions are rather close over the range 
shown in Fig.H. This means that the asymptotic behaviour of 
h{x), and more generally its behaviour for x > 10, typically, 
is still rather poorly determined. 



3.4 Moments of the counts in cells 

If one believes that the scaling function h{x) provides a rea- 
sonably good description of the counts-in-cells, it is readily 
seen from FigM that the measured counts in the simulation 
miss the very rare overdensities at a; > 10. In this region, 
actually there are counts that oscillate, due to (Bouchet et 
al.l992; Colombi et al.l995) whatever last cluster appears 
to be in the numerical data, reflecting their "cosmic vari- 
ance" . It turns out that the calculation of the moments Sp 
(m) badly reflects this missing information. The moments 
directly obtained from the measured counts are systemati- 
cally underestimated, the effect being more pronounced as 
the scale of the cells increases. The moments depend most 
on the counts-in-cells at a; ~ {p + LUs)xs, as can be seen from 
(h) and (H). With lOs — —I and Xs — 10 the calculation 
of S3, already, cannot be done without an extrapolation of 
the numerical data to the badly needed, but not measured 
large x (large N) tail. Appropriate methods have been pro- 
posed by the above authors , but the resulting values of Sp 
extracted from the data strongly reflect the precise extrap- 
olation procedure used by the various authors (Colombi et 
al.l996; Munshi et al.l999) . The differences get rapidly 



* To increase the available range, Colombi ot al (1996) correct 
for the oscillations (typically a factor of 3 from the mean, some- 
times an order of magnitude) due to this "cosmic variance" by 
means of an an educated guess, fitting a smooth curve of pre- 
determined shape to the data. To estimate the uncertainty, these 
authors determine by human judgement what can be considered 
as a reasonable fit and what cannot. This generates error bars 
given in the form of a factor /„ (i.e. Sp is within the range 

Sp,mean/fQ tO Sp^mean X /q) with typically fg = 1.1, 1.35, 1.65 

for p = 3,4,5, respectively, whereas the deviations from scaling 
show, under the same conditions, a scale-dependence by a factor 
/ = 1.2, 1.5, 1.85, respectively. Their interpretation is that, since 
/ is typically larger than fq by 10%, this " prooves a small but 
significant departure from the stable clustering predictions" in 



very large with increasing order p, from factors of 2 for p = 4 
to already an order of magnitude for p = 5. Thus extreme 
care must be taken in the interpretation of the Sp parame- 
ters extracted from the simulations. Rather than introducing 
sophisticated corrections to extract the large A'^ tail of the 
counts (with uncontrolled errors due to the badly needed 
correction for cosmic variance) it may be as accurate (inac- 
curate) to simply flt a function h{x) to the numerical data 
and use (0) to get the moments. 



4 MASS FUNCTIONS 

For astrophysical applications, one is usually more inter- 
ested in the mass functions (the number density of objects 
deflned for instance by a given density threshold) than in 
the counts in cells statistics. Moreover, one may wish to 
study a wide variety of objects, from low-density Lyman- 
Q clouds (which may even be underdense!) to clusters and 
massive very dense galaxies (see for instance Valageas & 
Schaeffer 1999 and Valageas et al.l999 for a detailed descrip- 
tion of Lyman-Qf clouds and galaxies in this framework). Of 
course, these very different objects are deflned by speciflc 
constraints so that they cannot be described solely by the 
mass function of "just-collapsed" halos with the traditional 
overdensity A — 177. 

We have determined numerically the multiplicity of var- 
ious kinds of idealized objects and compared them to the- 
oretical predictions. To this purpose, we use two different 
methods. The fist is the "spherical overdensity" algorithm 
(Lacey & Cole 1994) which ranks particles in order of de- 
creasing density and counts halos defined by the density 
contrast A by looking down this list (this introduces some 
double counting since part of a halo may be counted again 
at a lower rank: to avoid this, when a new halo is recognized 
its particles are removed from the list). The second method 
is based on a friend-of-friend algorithm with a linking length 
b which we relate to A through: 



(1 + A) = 178 (^) 



(22) 



4.1 The Press-Schechter approximation for 
just-collapsed halos 

Most studies have focused on the mass function of "just- 
collapsed" (or "just-virialized" ) halos which are defined by 
a density contrast in the actual non-linear density field 
A ~ 177 (for Q = 1). Indeed, this is the quantity which 
is considered by the popular Press-Schechter prescription 
(Press & Schechter 1974). The idea of this formulation is 
to recognize in the early universe, when density fluctuations 

the form of a scale-dependence of the coefficients Sp . Whereas we 
think that such a procedure (see Colombi et al.l994) to extract 
information on the large N behaviour, although somewhat indi- 
rect, is reasonable in the absence of any better way, it is not free 
from systematic biases and the associated correction may not be 
achieved at the accuracy needed to back the claim of the former 
authors. Indeed, a subsequent work (Munshi et al.l999), using a 
different computation, shows that the numerical data for large x 
are at least consistent with the assumption of no scale-dependence 
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are still gaussian and described by the linear theory, which 
overdensities will eventually collapse and form virialized ob- 
jects. Thus, one associates to any linear density contrast 5l 
obtained by the linear theory {Sl oc a for a given fluctuation 
of fixed mass) a non-linear density contrast A provided by 
the spherical model. Next one assumes in addition that col- 
lapsed objects virialize at a radius equal for instance to one 
half their turn-around radius at the time when the spheri- 
cal dynamics predicts a singularity. Then one obtains that 
when Sl reaches the threshold Sc — 1.69 a virialized halo 
with a density contrast A ~ 177 has formed. Finally, one 
identifies the mass fraction Fa(> M) within such objects 
of mass larger than M with the fraction of matter which is 
above the threshold 5^ at scale M in the linear universe: 



FA{>M) = FLi>5,;M) = 



dS 



-«^/(2ct^) 



s^ V2^a{M) 



(23) 



Then, the mass-derivative of the previous quantity provides 
the mass function of virialized halos: 



KM) 



dM 
M 



dM 



V2 



dv 



with 



a{M) 



(24) 



(25) 



Here we defined ^{M)dM/M as the fraction of matter which 
is enclosed in objects of mass M to AI + dM. Note that in 
(M) we have arbitrarily multiplied the mass function by the 
usual factor 2 so as to get the proper normalization to unity. 
We must however emphasize that this normalization correc- 
tion is not justified in the present case. Indeed, although this 
multiplicative factor 2 was recovered rigorously by Bond et 
al.(1991) using the excursion sets formalism for a top-hat 
in fc, taking into account the cloud-in-cloud problem, this 
result does not apply to more realistic filters like the top- 
hat in R used here. In fact, as shown in VS (and noticed 
by Peacock & Heavens 1990 through numerical results) this 
correction factor goes to unity at large masses. So, this fac- 
tor is not constant (and equal to 2) but (/-dependent. This 
simple normalization procedure (M) , nevertheless, leads to a 
scaling-law in the parameter v. the mass function ^{v)dv/u 
does not depend any longer on the initial power-spectrum. 
We present in FigJT^ the comparison of our numerical results 
with the PS prescription. We can check in the figure that 
the mass functions obtained by both methods are consistent. 
Note that the numerical points correspond to averages over 
different output times (weighted by the number of halos) 
of the mass functions realized in the simulations. Although 
there is some scatter we checked that all curves superpose 
on the mean mass function shown in Fig.R We can note that 
our results are consistent with the scaling in v obtained from 
the PS approach. 

Thus, we can see as was already checked by many au- 
thors (e.g. Efstathiou et al.l988; Kauffmann & White 1993; 
Lacey & Cole 1994) that the Press-Schechter mass func- 
tion (shown by the dashed line) agrees reasonably well with 
the numerical results. Indeed, although there are some small 
but significant discrepancies (especially at the low mass end 
where the slope seems to tend towards the non-linear predic- 
tion) . This could be expected in view of the simplicity of this 
approach but the agreement is still rather good for such a 




log(z^) 



log(i') 





log(jy) 

Figure 7. The mass function of "just-collapsed" halos defined by 
the density threshold A = 177, displayed in terms of the variable 
u. The crosses show the results of a "spherical overdensity" algo- 
rithm while the squares correspond to a "friends-of-friends" pro- 
cedure. The numerical results are averages over different output 
times from the simulation. The short-dashed curve is the standard 
PS prescription (note it does not depend on n). The dot-dashed 
line corresponds to ([2l|), referred to as PSs in the text, which is 
the same as the PS result in the case that the stable clustering 
regime is reached (and it can be seen that this is not fully the 
case here). The lower solid line is the scaling function h{x) ob- 
tained from counts-in-cells while the upper solid line is 8/7 h{x). 
The scaling model predicts that in the stable clustering regime 
the counts are in-between these two curves at large v, and below 
the upper one but with generically the same slope at small v. 
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simple model. A drawback of the Press-Schechter approach 
is that it is clearly limited to the mass function of "just- 
collapsed" objects and does not provide by itself predictions 
for mass functions defined at other density contrasts. 

4.2 The scaling model for the halo multiplicity 

Predictions for the mass function of objects defined by an 
arbitrary density contrast A in the actual non-linear density 
field can be obtained within the framework of the scaling 
model, as shown in detail by VS. These predictions are valid 
for any value of (H- A), whatever large or small (even below 
unity!), provided the conditions (H) are met. As was the case 
for the counts in cells the mass function exhibits a scaling- 
law in the parameter x: 

t^{M)— = x^ H{x) — (26) 

where we still have a; = (1 -I- A)/^(_R). Note that here R is 
the actual size of the halo (as opposed to the "mass scale" 
Rl, the radius an object of the same mass would have in a 
uniform density universe, which enters for instance in the PS 
mass function). The scaling function H{x) should be related 
to h{x) introduced for the counts in cells. Indeed, if one 
identifies the mass within objects (defined by the density 
threshold A) larger than R (i.e. more massive than M = 
(1 + A)p(47r/3)-R'^) with the fraction of matter enclosed in 
cells of scale R with a density contrast larger than A, one 
obtains a simple estimate Hceii of H{x) which reads 



Hceu{x) = h{x). 



(27) 



Even if this is probably an oversimplified approach, it 
has been steadily argued (Balian and Schaeffer 1989a ; 
Bernardeau and Schaeffer 1991; VS) at an increasing level 
of sophistication that the functions H{x) and h{x) can be 
expected to be close and that for practical purposes h{x) 
can be used as a good approximation to H{x). Also, as can 
be seen in VS, this procedure is found to be very close to 
the global derivative used in the PS prescription (E4|). In the 
scaling model, however, it is applied to the actual non-linear 
density field rather than to the original Gaussian field, and 
thus is correctly normalised. 

On the other hand, if one defines directly objects of scale 
R to R + dR by the constraints that the density contrast is 
larger than A on scale R but smaller than A on scale R+dR 
(thus one tries to recognize individually each halo in order to 
handle the cloud-in-cloud problem) one gets a new scaling 
function Hy^{x) which can be shown (VS) to satisfy (for 
constant A): 



Wx: Hy<{x)<-h{x) 

7 



x^ Xs : h(x) < Hy^{x) < — h(x) 

7 



(28) 



with generically (unless the leading order in the expansion 
accidentally vanishes) the same asymptotic behaviour at 
small x: 



a;< 1 



Hx{x) oc h{x) oc 



(29) 



Moreover, for a tree-model where the many-body correlation 
functions can be expressed as products of the two-point cor- 
relation function one obtains iif><(x) — h{x) for x ^ x^- 



Note that for a gaussian random field the same procedure 
(that is to define the mass function by requiring a given 
overdensity in the linear regime at a given scale R and a 
lower one at a slightly larger scale, in order to recognize in- 
dividually each object), would lead to a mass function which 
diverges at small masses (see VS), meaning divergent total 
mass due to the severe overcounting. Even for a fixed density 
contrast A, the same mass is counted many times at smaller 
and smaller scales because of numerous density inversions. 
Since the fraction of matter (compared to unity) enclosed 
in objects defined by a density threshold is a direct measure 
of the severity of the "cloud-in-cloud" problem, this shows 
the latter is indeed severe for Gaussian fiuctuations. On the 
other hand, in the non-linear case when one directly counts 
the overdensities in the non-linear density field, the normal- 
ization is necessarily lower than 3/7 for the case described 
above from (ESI) and can be expected to be close to unity. 
This shows, as argued by VS, that in the non-linear regime, 
the overdensities being well-defined, there is no longer any 
serious cloud-in-cloud problem. The latter still manifests it- 
self in the sense that the use of ff>< induces some slight 
overcounting, typically of the order of 20% (see VS). This is 
the same overcounting as the one found in computations for 
the spherical density algorithm (which has been removed in 
this case by attributing all the mass to the larger object), 
and indeed of similar magnitude, as we checked in the nu- 
merical simulations. The true constraint, rather than the one 
on the edge of the halos used to define H^ < , which avoids 
this overcounting, is the condition that the density contrast 
is smaller than A for all scales larger than R+dR. This would 
completely solve the cloud-in-cloud problem. However, the 
previous formulation leading to Hy^{x) should already pro- 
vide a satisfactory tool. Moreover, it is close to the algorithm 
actually used in numerical simulations. 

The solid curves in Fig.0 represent the functions h{x) 
and 3/7 h{x). From the results described above, we expect 
the mass function obtained from the numerical simulations 
to lie below the upper curve for all v and between these two 
curves for large u (which corresponds to large x) beyond the 
exponential cutoff. We can see that this is indeed approxi- 
mately the case for all power-spectra. However, it appears 
that the numerical mass function is significantly different 
from h{x), especially at the low mass end. This could be 
expected since it does not seem possible to derive a lower 
bound for H{x) in this domain in the general case. Never- 
theless, the slopes of both scaling functions appear to be 
the same H{x) oc h{x) for x ^ Xs- Note that for large v 
and X (x > A/10) the mass function should not behave 
as (bq) since this domain corresponds to large scales which 
are no longer in the non-linear regime, so the correlations 
functions do not obey the scaling (0). The dot-dashed lines 
correspond to the "PSs" mass function ( |26[ ) given by the 
scaling approach with H(x) taken equal to the "spherical 
model" scaling function (El|) . Of course at small scales it su- 
perposes onto the standard PS mass function, as explained 
in Sect. 3.3, since in this regime it is equivalent to the PS pre- 



diction. At large scales it differs from the latter because one 
leaves the highly non-linear, stable-clustering, regime and 
the two-point correlation function is no longer given by the 
asymptotic behaviour 1{R) = [10/(3a) o{RL)f , see (|l|). 
Note that the PS approach works somewhat better for the 
mass function than it did for the counts in cells statistics. 



© 0000 RAS, MNRAS 000, 000-000 



12 P. Valageas, C. Lacey and R. Schaeffer 



and that, for the contrast A = 177 relevant to the present 
calculations, at large masses the difference with respect to 
the stable clustering PSs result somewhat helps. What helps 
also to bring the PS result in agreement with the numerical 
data is the factor of 2 by which this mass function has been 
multiplied, as is traditional, despite numerical calculations 
(Peacock & Heavens 1990) as well as theoretical considera- 
tions (VS) show this factor should be close to unity at large 
X (and is expected to be larger at small x). 

As compared to the standard PS approach, the inter- 
est of the formulation (|26| ) for the mass function of "just- 
virialized" objects is that it makes a connection with an- 
other characteristics of the density field: the counts in cells 
statistics described earlier. 
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4.3 Different constant density contrasts 

The main advantage advantage of the scaling model is that 
it can deal with more general mass functions. We shall first 
consider the case of mass functions which are still defined 
by a constant density contrast but where A can now take 
any value. In the scaling approach which led to l)^) the 
quantity (1-1- A) is only used to define objects but no specific 
value is singled out as opposed to the PS mass function 
which explicitly deals with "just-collapsed" objects (so that 
it makes sense for only one class of halos characterized by 
a non-linear density contrast A ~ 177). Hence the relation 
(M) should hold for any value of (1 -I- A), in the domain 
where the scaling (0) holds. 



4-3.1 Spherical overdensity agorithm 

We show in Fig.g the mass functions we obtain from the nu- 
merical simulation for various values of the density contrast 
for the "spherical overdensity" algorithm. We display the 
quantity n{M) dlnM/dlnx as a function of a; = {1 + A)/^{R) 
(as explained previously ^(-R) was obtained from the same 
simulation through counts in cells, see Fig.H). If the scaling 
(M) holds all curves should superpose (in the relevant range 
of validity of the scaling laws) onto the function x^ H(x) 
(we only show the range s < (1 -I- A) which corresponds to 
^ > 1 but the range of validity could be smaller than this: 
^ S> 1). The points correspond again to an average over 
different output times. 

We can see that the available counts are pushed much 
further towards the larger values of x as compared to Fig.W 
This is because we check the neighbourhood of particles to 
get the halos and thus we are guaranteed to be in a re- 
gion where there really is some mass. To get the statistics of 
the counts-in-cells, on the other hand, implies to check ran- 
domly over the whole volume, so that the statistics of the 
less dense but more extended regions are r elati vely favoured. 

this amounts 



As will be discussed in more detail in Sect, 4.4, 
in some sense to measure xH{x) rather than H{x), with a 
better statistics at the large values of x (note, however, that 
whereas x^H{x) has a maximum, both the above functions 
are decreasing for increasing x, still favoring the small val- 
ues of x as compared to the large ones) . A similar effect was 
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described in Sect. 3.2 for the measure of ^ as opposed to ^. 
As expected, all curves lie between h{x) and 3/7 h(x) 
beyond the exponential cutoff and are located below 
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log(x) 

Figure 8. The mass functions of halos defined by various den- 
sity thresholds A obtained from a "spherical overdensity" algo- 
rithm. Different symbols correspond to different values of (1 + A). 
The lower solid curve is the scaling function h(x) measured from 
counts in cells while the upper solid curve is 3/7 hi^x). The non- 
linear model predicts that the counts are in-between these two 
curves at large x, and below the upper one but generically with 
the same slope for small x. 



3/7 hix) for all values of x. Moreover, the slope at small 
X is in most cases close to that for h(x), although for n = 0, 
where the statistics are the best, a distinct difference is ap- 
parent: whereas h{x) has a slope (Tab. 2) oj = 0.45, here the 
slope is rather u) = 0.6. 

Moreover, although all curves superpose beyond the 
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cutoff as predicted, there appears to be a shift in the power- 
law domain with (1 + A): larger density contrasts lead to 
a smaller amplitude in this part of the mass function. This 
drift appears at a; < 3 for (1 + A) = 5000 and x < 1 for 
(1 + A) = 1000 and has disappeared for (1 + A) = 178 (in- 
deed the curves obtained for 10 < (1 + A) < 178 superpose 
very well). 

These deviations will be discussed in the following sec- 
tion. 



4-3.2 Friend-of-friend algorithm 

We show in Fig.y the mass functions we obtain from the 
"friend-of-friend" algorithm. We can check that they are 
consistent with the results from the "spherical overdensity" 
algorithm. For large A we find as previously that the curves 
superpose beyond the exponential cutoff but that there is a 
drift in the power-law domain. 

The advantage of the "friends-of-friends" algorithm is 
that we can also consider negative density contrasts, provid- 
ing for much larger a range for the possible variations of A. 
In such a case we count one large halo which extends over 
the whole box (the clusters linked by the filaments create 
a percolating network) and also small underdense objects 
which appear as density peaks within voids and are sepa- 
rated from the main cluster by their very low density sur- 
roundings. Thus, the mass function defined in this way only 
makes sense for very small masses (small parameter x). De- 
spite the low density, this is still in the strongly non-linear 
regime since a; ^ (1 + A) so that ^ S> 1. 

For n = 0, we can see that all curves fall on top of 
each other, down to (1 + A) = 0.5 which is a little low at 
logx < —0.5. Similarly, for n = — 1, the results exhibit the 
scaling theoretically implied by the stable clustering, except 
for (1 + A) — 0.5 at the larger values of a; (log a; > —1.5). For 
n = — 2, also, the curves with (1 -I- A) = 0.5 and (1 + A) = 1 
do not fall on the scaling curve built by the result at larger 
density contrasts. As shown in Fig.hd, the deviations found 
for the lower values are seen to be closely associated to the 
degradation of the statistics (the counts are low at the larger 
a;), that occurs for all three values of n we have studied 
precisely at the same values of (1 -I- A) and log a: where the 
deviations are noticed. Moreover, the deviation at large x 
of the mass functions measured for low density contrasts 
(A < 0) from the scaling predictions corresponds to the fact 
that they only make sense for small x since we must have 
a; <C (1 -I- A) so as to avoid entering the linear regime where 
there are no dense spots within underdense regions like the 
ones we look for here. 

So, the non-linear scaling prediction is seen to hold 
down to values of (1 + A) that are unity or even smaller. This 
is nearly three orders of magnitude below the standard den- 
sity contrast of 178 usually studied. It has also been tested 
down to log (a:) = —2.5, a real improvement as compared 
to existing work; values below log (a:) = — 1 had never been 
considered previously. 

Note that the slopes at small x of the mass function 
are close to the ones predicted by the theory. In the n — 
case, for instance, u ~ 0.5 , quite close to the slope to — 
0.45 of h{x). The offset seen in Fig.y with the "spherical 
overdensity" algorithm does not exist in this case. It should 
be noted that the n = case with (1 -)- A) — 0.5 clearly 
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Figure 9. The mass functions of hales defined by various density 
thresholds A obtained from a "friend-of-friend" algorithm. Diff'er- 
ent symbols correspond to different values of (1 + A). The data 
shown correspond to averages over several simulations. The lower 
sohd curve is the scahng function h{x) measured from counts in 
cells while the upper solid curve is 3/7 h{x). The non-linear model 
predicts that the counts are in-between these two curves at large 
X, and below the upper one but generically with the same slope 
for small x. 
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Figure 10. The mass functions obtained from a "friend-of-friend" 
algoritlini for (1 + A) = 0.5 (witli n = and n = —1) and 
(1 + A) = 1 (witli n = —2), as in Fig.H. Different symbols cor- 
respond to different ranges for tfie comoving size R of the halos. 
The various points located at the same value of x for a given 
range of R correspond to different times (and slightly different 
R within the allowed range). They should superpose due to the 
self-similarity induced by the power-law initial conditions. 

shows (FigJ9|) that the non-linear power-law prediction (and 
perhaps more important, the prediction of existence !) for 
objects that are much denser than their close environment 
is well verified even in this quite extreme case of very low 
density, definitely below the average. 

Towards the higher densities, we have run calculations 



up to (1 + A) — 5000, a factor of 30 above the usual den- 
sities. As for the spherical density algorithm, scaling is well 
verified. Similarly, at the smaller values of x, steadily in- 
creasing deviations appear when A is increased, although 
to a lesser extent. For n = —1 and log(a::) — 0, the curve 
for (1 + A) = 5000 is below the one for (1 + A) = 1000 
(that at fixed x corresponds to correlation functions 5 times 
smaller, hence to scales 5^'^'^ = 2.9 times larger at fixed 
time) by a factor 1.35 and below the (1 + A) — 178 curve by 
a factor 1.7. The relevant scales are the smallest reachable 
in the simulation, typically R = 0.5 Mpc at log(a;) = and 
R < 0.3 Mpc at log(a;) = -0.4, for (1 + A) = 5000. This is to 
be compared with the force softening parameter which is 0.2 
Mpc in all the calculations presented here. We have checked 
the sensitivity of the calculations to variations of this soft- 
ening parameter. Decreasing the latter by a factor of 2 so 
as to increase the effect of collisions, lowers the counts by a 
factor of 1.25 at the scales considered above (a change quite 
comparable to the deviations from scaling that were found) , 
but has no effect at the larger values of x (the larger scales). 
For the calculation of the correlation function, at ^ > 1000, 
similarly, deviations of the computation appeared at small 
scales. So, there is a limit towards the large densities and 
the smaller values of x beyond which the numerical tests 
are no longer conclusive, since binary collisions due to the 
discreteness of the points in the simulation start to play a 
non-negligible role. 

Although the previous discussion leads to attribute the 
lack of objects below some value of x, for a given (large) den- 
sity contrast, to the fact that the smallest available scale in 
the simulation has been reached, this is an important point 
to be discussed. Note that the theoretical prediction (see 
VS) for H{x) extends to all values of x (including the lower 
ones for which there are no numerical data) and that, what- 
ever the density contrast, the area under the scaling curve is 
unity (it is actually slightly above for 7y>< (a;), but is unity if 
corrections to avoid double counting are made) . This means 
that however large the required density threshold is, it de- 
fines a way to distribute all the mass among objects having 
this density contrast. The theoretical prediction thus calls 
for very strong suhclustenng at all density levels. On the 
other hand, the simulations do not provide enough power at 
small X for large overdensities: there is barely more than 50% 
of the mass in objects of contrast 5000, for instance, showing 
an increasing lack of high density objects. Again, this ap- 
parent difference is worth to be thoroughly investigated in 
future more extended simulations, to see whether subclus- 
tering up to extremely high overdensities is present down 
to very small scales or not. However, as shown in Valageas 
(1999), note that the reasonable agreement of the numerical 
results with our scaling model already shows that there is a 
large amount of subclustering. It has been steadily argued 
(Balian & Schaeffer 1989a,b) that the existence and the be- 
haviour of these substructures are the base of the scaling 
model and govern the properties of the density field. 

We have tested the scaling model for the multiplicity 
function for vastly different density contrasts, spanning four 
orders of magnitude, from (1 + A) = 0.5 to (1 + A) = 5000, 
finding non-linear objects at every density. We have shown 
that up to the small-scale limit where collision effects be- 
gin to play a role and to the low-density limit where devia- 
tions from self-similarity start to appear the scaling model 
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provides a reasonable approximation to the mass functions 
obtained in the simulations. Also, we can note that the scal- 
ing model gives surprisingly good results in domains which 
are beyond its expected range of validity. For instance, the 
curves for (1 + A) = 10 and (1 + A) = 178 agree in the 
numerical results down to ^ ~ 1. Hence, although the model 
may not be perfect it is still quite powerful as it allows one 
to get a good estimate of the mass functions of very dif- 
ferent objects. Moreover, it clarifies an interesting link with 
another statistical analysis of the density field, namely the 
counts in cells. 



4.4 Constant radius 

In practice, one may also be interested in objects which are 
no longer defined by a constant density threshold but by a 
condition of the form (1-1- A) oc R~^^ . The predictions of the 
scaling model for such a case were studied in VS. It hap- 
pens that for astrophysical objects like galaxies or Lyman-a 
clouds the constraints which define these mass condensations 
can be approximated in some range by the requirement that 
their radius is constant and equal to a given scale which may 
be associated in the first case to a cooling radius and in the 
second to a Jeans length (see Valageas & Schaefi'er 1999 and 
Valageas et al.l999 for more details). This corresponds to 
the limit /3 — > oo. Then, one has to add the supplementary 
constraint that the density profile is locally decreasing (since 
one usually wants to consider peaks and not valleys). In this 
specific case one obtains the relations (see VS): 



Hceu(x) = h{x) 



(30) 



for the simple model. More realistically, one can define (see 
VS) a scaling function H^^ix) for the mass distribution of 
objects of fixed radius and make the statistics of the overden- 
sity (1 -I- A) -hence of the mass since we have M oc (1 -I- A)- 
by requiring for such a halo to have a density contrast larger 
than A at scale R and smaller than A at a slightly larger 
scale. Thus one considers individually each object, defined 
by the fixed radius R and its density contrast A (which gives 
its mass), making sure that it is surrounded by regions of 
lower density. Then, one obtains the bounds 



Va; : ii'><(a;) < h{x) 
x^^Xs : ii"><(a;) > 7/3 /i(a;) 



(31) 



Obviously this procedure is quite similar to the counts in 
cells so we can expect the scaling function H{x) to be very 
close to h{x). 

We show in Fig.hll the mass functions obtained in this 
way from the numerical simulations. We simply use a mod- 
ified version of the "spherical overdensity algorithm" , look- 
ing at particles in order of decreasing density and defining 
halos as objects of constant size R around density peaks. 
The scaling with the variable x is well verified, for all values 
of the spectral index n, with a result that is very close to 
the prediction dSQ). There seems, nevertheless, to be more 
counts at large :r(large masses) : the first of the bounds ( pjl[ ) 
of the "improved" estimate -H'><(x) (in the sense it is closer 
to the spherical overdensity procedure) seems to be some- 
what violated. Note however that the curves obtained for 
different radii (and times) superpose fairly well which shows 
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Figure 11. The mass functions of halos defined by various co- 
moving radii /?, obtained from a "modified spherical overdensity" 
algorithm (see main text). Different symbols correspond to differ- 
ent values of R. Note that for each radius we display the results 
obtained at several times which correspond to various values of 
5. The upper solid curve is the scaling function h(x) measured 
from counts in cells while the lower solid curve is 7/3 h(x). The 
non-linear model predicts that the counts are in-between these 
two curves at large re, and below the upper one but generically 
with the same slope for small x. 
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the predicted scaling in x is verified (within the numeri- 
cal scatter). In fact, this small deviation is probably due to 
the definition of _H"><(i) itself and not to a failure of the 
scaling-laws (12). Indeed, this estimate of the mass function 
is derived from specific counts-in-cells (with the characteris- 
tic that these cells consist of an internal sphere surrounded 
by a small corona) so that one can miss the most massive 
high density halos by looking at a sphere which is not ex- 
actly centered onto the underlying halo. This decreases the 
mass enclosed in this sphere, which is maximum for a cell 
correctly centered. On the other hand, in the simulation one 
directly draws the spheres from the highest density peaks, so 
that one always counts the largest amount of matter which 
can be attached to a given halo (if the object is approxi- 
mately spherically symmetric). This latter procedure (look- 
ing at the peaks) is not included in _ff><(2;) and, because 
of this small "mis-placement" of the spheres, it slightly un- 
derestimates the number of very high density and massive 
halos. We note that a similar, but very small, deviation ma- 
be seen for the mass functions shown in FigJS and Fig 
for the same reason. We also note that the mass function 
shown in FigJllI allows one to probe deeper into the expo- 
nential cutoff of H(x), hence of h{x), as compared to the 
counts in cells shown in Fig.W. This could also be expected 
from the fact that in this procedure "cells" are drawn di- 
rectly around highest density peaks, which are thus well ac- 
counted for, while for the statistics of P{N) density peaks 
can be divided between several cells. 



5 CONCLUSION 

In this article we have presented a comparison of the results 
obtained from numerical simulations to the analytical pre- 
dictions implied by the Press- Schechter (Press & Schechter 
1974) prescription and the scaling model (Balian & Schaeffer 
1989a ; VS). 

We have first checked that the two-point correlation 
functions follow the behaviour predicted by the stable- 
clustering ansatz, as shown by other studies. However, al- 
though there is a qualitative agreement between various 
works the exact value of the amplitude of the two-point 
correlation function in the highly non-linear regime still re- 
mains poorly determined as there remains some discrepancy 
between different studies. Next we have shown that the scal- 
ing model provides a good description of the counts in cells 
statistics. The characteristic scaling functions h[x) we ob- 
tain over the unprecedented range —2.5 < log x < 1 agree 
reasonably well with previous estimates. Nevertheless, the 
asymptotic behaviour of the exponential cutofi' is not very 
well constrained for n — —2 due to the limited range of 
scales available in the simulations. 

Then we considered the mass function of "just- 
collapsed" objects, which is the quantity most authors have 
focussed on. As was already noticed by many studies, we find 
that the PS approximation works reasonably well for density 
contrasts of ~ 200 (although there are some discrepancies) . 
On the other hand, the numerical results are also consistent 
with the predictions of the scaling model (although they do 
not provide the exact value of the mass function) . Next we 
have studied more general mass functions defined by vari- 
ous density contrasts that vary by four orders of magnitude 



(including negative density thresholds!). The scaling model 
(which is the only currently available analytical tool to han- 
dle these quantities) was shown to provide a reasonable es- 
timate of these mass functions. 

Finally, we considered the limiting case of objects de- 
fined by a constant radius constraint (which arises in an as- 
trophysical context from cooling conditions). We have shown 
that the scaling model also works very well for such mass 
functions. Some small deviations appear in certain regimes. 
We make specific suggestions on which biases may appear in 
the simulations and show they are of the right order of mag- 
nitude. Although these biases likely explain the observed 
scatter, it is quite clear that we cannot check the scaling 
prediction better than this scatter. It remains to be verified 
in more powerful simulations whether (at least part of) these 
deviations can be attributed to some violation of the scaling 
predicted by theory, or if these deviations indeed disappear 
and thus simply refiect the fact that we have pushed the 
present simulation as far as possible, and reached its limits. 

Thus, the scaling model provides a reasonable descrip- 
tion of the density field in the non-linear regime (although 
there are some deviations in certain regimes they remain 
reasonably small for practical purposes). Moreover, it allows 
one to link two different properties of the latter: the counts 
in cells statistics and the mass functions. In addition, it can 
be used to obtain many different mass functions (in addition 
to the usual "just-collapsed" halos) which is of great interest 
for practical purposes when one intends to model various ob- 
jects like Lyman-a clouds or galaxies which clearly cannot 
be defined by the sole constraint A = 177. Obviously, our 
study should be extended to more realistic power-spectra 
which are no longer power-laws (e.g. CDM) and to different 
cosmologies (low-density universes). 
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